A phase-field model of Hele-Shaw flows in the high viscosity contrast regime 
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A one-sided phase-field model is proposed to study the dynamics of unstable interfaces of Hele- 
Shaw flows in the high viscosity contrast regime. The corresponding macroscopic equations are 
obtained by means of an asymptotic expansion from the phase-field model. Numerical integrations 
of the phase-field model in a rectangular Hele-Shaw cell reproduce finger competition with the final 
evolution to a steady state finger. 
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The characterization of the dynamics of morphologi- 
cally unstable interfaces is one of major problems of non- 
equilibrium phenomenology [1]. Some relevant examples 
of interfaces that grow out of equilibrium are dendritic 
growth, directional solidification, flow in porous media, 
electro-deposition, bacterial colony growth and two-fluid 
flow in a Hele-Shaw cell. The latter example is also called 
the Saffman- Taylor problem and has played a central role 
in this field, both because of its relative simplicity and 
because of its potential importance in oil recovery. It has 
been widely studied both experimentally and theoreti- 
cally. 

Even if the Saffman- Taylor problem is mathematically 
simple in relation to other problems, it has a moving 
boundary condition which makes it a free-boundary prob- 
lem. The corresponding equations have been solved ana- 
lytically for very short times by means of a linear stability 
analysis and for the steady state finger shape by means of 
conformal mapping techniques [2,3]. Some analytical re- 
sults have also been obtained for the dynamics of interme- 
diate time [4]. Numerically there are several techniques, 
most of them involving integral boundary methods [5-7] 

The so-called phase-field models have been introduced 
within the context of solidification to study the dynam- 
ics from the linear regime to the long time behaviour [8] . 
These models are based on the introduction of a meso- 
scopic equation for an order parameter (the phase- field) . 
This equation is coupled to other physical fields (such as 
a thermal field). The advantage of this method is that 
one does not have to explicitly trace the interface. It is 
a field model for all values of the order parameter that 
varies continuously from one phase to the other. One 
has to identify the locus of points with a given value of 
the order parameter, which is arbitrarily chosen to be 
the interface. The use of a mesoscopic model, for which 
the interface has a small width e, is justified as long as 
in the sharp interface limit, e ^ 0, the correct macro- 



scopic equations are recovered. Recently, the concept of 
phase-field models has been used in a broader sense to 
include any model which contains continuous fields that 
are introduced to describe phases separated by diffuse 
interfaces. Phase-field models have been used in a wide 
range of problems such as viscous fingering, roughening, 
vesicles, pinchoff and reconnection in a Hele-Shaw cell 
and intracellular dynamics [9-14]. 

In general, the phase- field models have been considered 
for symmetric situations where the characteristic param- 
eters (such as the thermal diffusivity) are identical in 
both phases. This gives rise to the so-called two-sided 
symmetric models. Very recently. Karma [15] has pro- 
posed a phase-field model of the one-sided type (with 
zero diffusion in one phase) to simulate quantitatively 
micro-structural pattern formation of alloy solidification. 
For the viscous fingering problem with arbitrary viscosity 
contrast, a phase-field model has been introduced in Refs. 
[9] . Such phase-field model, which is a two-sided model, is 
useful to describe the problem of viscous fingering except 
in the high viscosity contrast regime. This regime is ex- 
perimentally relevant since typically the pushing fluid is 
air or other fluid of negligible viscosity. For such a regime, 
a proper model was laking and this is what we are pre- 
senting in this paper: a one-sided phase field model for 
the high viscosity contrast regime of the viscous fingering 
problem. 

Our model contains an equation for an order parame- 
ter. It is Model B of Ginzburg-Landau phenomenology 
[16]. Instead of the coupling of the order parameter to 
a physical field through a second equation, we include a 
boundary condition such that the interface becomes un- 
stable. This is done by means of a ramp that creates a 
flux from the boundary. To consider a one-sided model, 
we only need to neglect changes in the order parameter 
in one of the two phases. The model could also be rele- 
vant for dendritic growth at very small undercooling by 
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introducing anisotropy [17,18]. 

Our phase-field model has the advantage of being very 
simple to be implemented on a computer and contains 
a complete description of all the nonlinear and nonlocal 
properties of the macroscopic model. We show how the 
macroscopic equations of the problem are obtained from 
the phase- field model in the sharp interface limit. This 
is done by means of the matched asymptotic expansion 
method. We then present numerical solutions showing 
how our phase-field model reproduces the main features 
of the viscous fingering problem such as the dynamic 
competition of modes and the formation of a steady state 
finger. This makes the model an attractive tool to be 
used to study problems that would not be easily feasi- 
ble with traditional methods such as the propagation of 
viscous fingering in the presence of quenched disorder. 



THE MODEL 

The Viscous Fingering Problem 

In the Saffman Taylor problem both fluids are gov- 
erned by Darcy's law, which relates the fluid velocity to 
the pressure gradient [2]. When the low viscosity fluid 
displaces the high viscosity fluid, the interface between 
both fluids is unstable. When the pushing fluid is con- 
sidered to have zero viscosity, Darcy's law states that the 
pressure on the pushing fluid is constant and all that re- 
mains to be solved are the equations for the viscous fluid 
subject to the proper boundary conditions at the fluid- 
fluid interface. This is called the high viscosity contrast 
regime and in this regime the equations for the displaced 
viscous fluid are 



V^p = 0, 



Vn = -KVp ■ n, 



(1) 



(2) 



the equations for our phase-field model and show how 
it reproduces the above equations in the sharp interface 
limit. 



Phcise-fleld model 

Our phase-field model contains a time-dependent 
Ginzburg-Landau equation for a conserved order param- 
eter and includes a boundary condition that makes the 
interface unstable. The equation reads 



f = V.[M„V,- 



(4) 



The local order parameter (f> adopts the equilibrium val- 
ues = 1 (air-phase) and 4>eq = — 1 (viscous fluid- 
phase). At the interface, varies continuously from one 
phase to the other. The parameter Mq has a constant 
value in each phase and is zero in air 



Mo = 



M, 



if 
if 



< 0, 
> 0. 



(5) 



VISCOUS 

fluid 



(3) 



Eq. (1) is the Laplace equation in the bulk, where p is the 
pressure of the viscous fluid. At the interface, there are 
two boundary conditions: the continuity equation, Eq. 
(2) and the Gibbs- Thomson condition, Eq. (3). v„ is the 
velocity normal to the interface. K is the permeability 
of the viscous fluid, K — where b is the separation 
between the plates and ji is the viscosity of the fluid that 
is being pushed. Ap is the pressure of the viscous fluid 
minus the constant pressure at the zero viscosity fluid, 
which without loss of generality can be taken equal to 
zero. K is the local curvature at the interface and 7 is the 
surface tension. These three equations also describe so- 
lidification in the quasi-static limit of small undercooling 
by introducing anisotropy. In what follows, we present 



FIG. 1. Scheme of the initial profile prepared with a ramp 
that will be maintained during the temporal evolution. 



The air phase can be pulled toward the viscous fiuid. 
An unstable interface is developed by maintaining a slope 
in the order parameter close to the interface, as the case 
shown in Fig. 1. This situation can be created by initially 
preparing the system with a profile of the form 



(l>{x,y) = 




if 2/ > 2/m, 

Urn) a Um - I < y < Vm, 

y < Vm - I, 



(6) 



and fixing the value 4>f = ~1 + al behind the interface 
up to a distance I throughout all the temporal evolution. 
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This slope represents the driving force of the system. The 
parameter a controls the slope and the fingers growth ve- 
locity. For convenience we will refer to the air as the plus 
phase and to the viscous fluid as the minus phase. 

THE SHARP INTERFACE LIMIT 

In this section we obtain the macroscopic equations for 

the viscous fingering problem in the high viscosity con- 
trast regime by means of an asymptotic expansion of the 
phase-field model in the sharp interface limit e ^ [19]. 

Our starting point is Eq. (4) used in the study of a 
conserved order parameter, 0. The chemical potential is 
given by 

fii(|))=^lB-e^V^cf> = -(^ + ^''-e^V^. (7) 

We divide the space into an outer and an inner region. 
We assume that <p = ±<peg + 0(e) far from the interface. 
e is considered to be a small parameter and we expand all 
the variables a(r, t) around the value e = in the outer 
region. We obtain 

a{r,t) = aa + eai + f^a2 + ... (8) 

For the interfacial region or inner region, we adapt our 
coordinate system using time dependent curvilinear co- 
ordinates. The interfacial points arc given by the curvi- 
linear coordinates u, which is the normal distance to the 
interface, and s, which is the arclength. Because the 
natural dimension in the inner region must be small , we 
introduce the variable w defined a.s w = ^. Thus, in the 
sharp interface limit, when e — »■ 0, the inner region goes 
from w ^ —oo to w +00. We use the corresponding 
inner fields A{w, s, t) in the inner region and the corre- 
sponding expansion is 

A{w,s,t) = AQ + eAi+e'^A2 + ... (9) 

When we take the limit of the sharp interface, e — > 0, 
the conditions for the field a and A , from the expansions 
to i-th order in e are 

lim Ai = lim a^, (10) 

w^ — oo w— »— 

lim dwAi+i = lim a„Oi. (11) 

w—*-—oo u— >— 

Due to the fact that m = in air, the matching condition 
is only imposed in the viscous phase. 

In the inner region, we introduce the order parameter 
^ such that <f){u{t), s, t) — 4>{r, t), therefore 

dtcj) = dt4> + dtudu4>- (12) 

We rescale time as t = et since we work in the quasi- 
static approximation, where the characteristic times for 



interface motion arc much larger than the characteristic 
times for the diffusion to take place. The local curva- 
ture K = — V^w is positive when a bump of the (/> > 
phase protrudes into the (f) < phase. Starting from Eq. 
(12), using the Laplacian operator in curvilinear coordi- 
nates = du^ — ndu + dg, and doing the corresponding 
variable changes, we have 

€dr<f> - -d^cf) = M{\^J^J^{(j)) - -a^M(0) + 9. Vl-/-)), (13) 

where we have dropped the tildes. The normal velocity 
V = —dtu is positive if the phase with a negative order 
parameter goes into the phase with a positive order pa- 
rameter. This variable is also expanded in powers of e. 

For the chemical potential n, the inner expansion in 
terms of (j) (to order e^) is given by 

/i((?i) = /io + eMi +e^A*2, (14) 

with 

Mo = M-Bo - dw4>o, 

Ml = Mbo*^! ^ ^^^1 + '•^■9w(f>o, (15) 
M2 = ^Mbo('/>i)^ + MBo<?^2 - dl(j)2 + Kd^cpi - 5f0o, 

where ~ Ms('/'o)- The prime represents the deriva- 
tive of (j) evaluated at (po. 

For the region far from the interface (outer region) , the 
length scale involved is much greater than e, so we can 
use a common time-independent coordinate system. In 
the viscous fluid region the dynamical equation for the 
order parameter is simply 

edr^ = MV^ni^), (16) 

where m(</') is given by Eq. (7). 

Inner region 

We now proceed to solve the equations for the inner re- 
gion Eqs. (13)- (15). We also use the matching conditions 
Eq. (10) and Eq. (11). Solutions that obey <?i(0) = and 
dw4'o{—oo) =0 are required. 

order For the inner region, the dynamical equa- 
tion to lowest order in e, (e~^) is taken from Eq. (13) 

a>o = 0. (17) 

Here we have taken into account the expansion for m- 
The previous expression has a solution /xq = mo -|- now. 
The requirement that mo must be finite for w — > — oo im- 
plies that no = 0. Finally, we consider mo = and then 
(j)eq = ±1. Therefore, iiq — in the inner region. 

order Taking the first-order terms from the 
dynamical equation in the inner region, Eq.(13), we have 
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(18) 



Outer region 



since /xq = 0. Integrating Eq. (18) in w we find 

-vo(l)o = Md^fii + m. (19) 

By evaluating between the limits ^ = —00 and w = 00 
Eq. (19), we obtain 

-2<j)eqVo = Mdyjiii{-oo), (20) 

where 2^eq is the order parameter change between the 
two phases and we have only the contribution of the vis- 
cous fluid phase on the right-hand side. Using Eq. (11) 
we have dyjiJ,i{—oo) = 9ti/xo(— 0) = 0. ^From Eq. (20), vq 
also vanishes and Eq. (19) gives rii = 0. 

By integrating Eq. (19) in w, we find that = 
/xi(— 00) is a constant. In order to obtain /ii(— 00), we 
use its expression from Eq. (15) and we multiply both 
sides of this expression by dyj(j)o and integrate in w 

/xi(-oo) j dwdyj(l)o = j dwdyjcpoin'so - 5^)'/>i 

+ dw{dy,<t>of. (21) 

The function dw(j)Q is known as the Goldstone mode and 
is related to the translational invariance of the interface. 
The equation for (/)o, written as a function of the rescaled 
variable w, is /ibo ~ dw't'o = 0. Differentiating with 
respect to w we obtain an equation for dyj(l)o, which is 
il^Bo ~ 9^)9w4'o = 0- So the Goldstone mode is a zero 
eigenvector of the linear operator /i'^^ ~<^w- doing in- 
tegration by parts, the first term on the right-hand side 
of Eq. (21) vanishes and we obtain 

Mi(-O) = -^K, (22) 

(Peg 

where 7 = ^/ dw{dw(f>o)'^ is the surface tension and we 
have used the matching condition for /ii(— 0) from Eq. 
(10). Taking into account the fact that at the interface 
p{—0) = (/)egA*i(— 0) we obtain Eq. (3). 

order e° In order to obtain the continuity equation, 
we need to go to the next order. The dynamical equation 
to order e° in the inner region is 

-vidw4>Q = M{du,^H2 - KdwiJ,! + ds^no). (23) 

Integrating Eq. (23) in the direction normal to w we 
find 

-2(/.e,wi = Md^H2{-^) = MduHi{-Q), (24) 

where we have used the matching condition Eq. (10), and 
the fact that 9u,/xi(— 00) = and that /xq = 0. Eq. (24) 
could be written as Eq. (2) in terms of the pressure at 
the interface of the viscous fluid, where K = Mj {2(j)1). 



order e ^ The dynamical equation in the outer region 
to lowest order is 

V^/xo = 0. (25) 

The boundary condition far from the interface is then 
/Uo = 0. We previously found that /Uq = for the inner 
solution at the interface. The only solution satisfying 
both conditions is = 0. It follows that 4)q = ct),,,j in the 
plus phase and 0o = ~</'eg in the minus phase. This was 
to be expected since the lowest order in the expansion 
corresponds to the solution of the fiat interface, 
order The dynamical equation for the order 

is 

V^/xi = 0. (26) 

At order e the order parameter and the chemical poten- 
tial are proportional and from Eq. (26) we obtain Eq. 
(!)• 



NUMERICAL RESULTS 

We have numerically integrated Eqs. (4)- (5) with e — 1 
and M = 1 on a rectangular lattice of vertical size 
Ly = 200 and mesh size Ao; = 1, with periodic boundary 
conditions in the x direction and refiecting boundary con- 
ditions in the y-direction. The system has been prepared 
with a horizontal interface containing some perturbations 
in order to be destabilized. The profile in the vertical di- 
rection is done by Eq. (6) with / = 10. As was mentioned 
before, during the evolution we maintain a slope by fix- 
ing the value (f>f = —1 + al behind the interface, up to 
a distance I measured from the tip of the most advanced 
finger. 



Finger competition 

Firstly, we arc interested in the generation and sub- 
sequent competition of fingers during the early stages of 
the evolution. A wide system of size = 128 has been 
considered and we have prepared an initial corrugated 
interface formed by the superposition of several modes 
of random amplitude. In Fig. 2 we show a typical evo- 
lution. It is seen that fingers develop from the random 
initial configuration. Some modes grow, some modes de- 
cay and finger competition begins. Both features have 
been observed in theoretical and experimental studies of 
the viscous fingering problem. The competition process 
continues until only one of the fingers survives. 
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(a) 



(b) 



(c) 



(d) 



FIG. 2. Finger development and competition for a = 0.04 
corresponding to early times: t = 100 (a), 500 (b), 1000 (c) 
and 2000 (d). 

In order to better visualize the competition process we 
have prepared a second initial condition consisting of two 
well-formed fingers, in which one of them is a bit larger 
than the other. In Fig. 3 we observe^ how the largest fin- 
ger grows at the expense of the other, which moves back- 
wards becoming smaller and eventually disappears. 
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(a) 



(b) 



(c) 



(d) 



FIG. 3. Finger competition process for two initially 

well-formod fingers with a = 0.04 and system width hx = 64. 
The patterns are separated by time intervals of 1000. 



Steady-state finger 

The width of the steady-state finger is expected to go 
to one half of the channel width as the velocity of the 
finger tip increases [3,21]. To better explore this sit- 
uation, we have considered a narrow channel of width 
Lx = 32, prepared with an initial condition that gives a 



single finger. We have analysed the temporal evolution 
of the finger for different tip velocities, corresponding to 
different values of the parameter a. In agreement with 
known results, higher velocities led to narrower fingers 
An example of the interface evolution is shown in Fig. 
4a. In Figs. 4b and 4c we compare the finger shapes 
obtained numerically for two different values of a with 
the theoretical shape for the Saffman- Taylor finger [2]: 



y = utip 



L.O- - A) 
27r 



In 



1 + cos 



7r(2a; - L^) 



(27) 



being A the ratio of the width of the finger to the width 
of the channel. To determine A from our numerical re- 
sults, we have evaluated the average width of the finger 
throughout the evolution, in a strip of thickness e = 4 
placed at distance 40 from the tip. For high enough tip 
velocities, our numerical results are in agreement with the 
Saffman- Taylor solution since they correspond to values 
of A closed to 1/2 (Fig. 4c), where surface tension ef- 
fects arc negligible. Also, the expected deviation from 
the Saffman- Taylor solution is observed for wider fingers 
in qualitative agreement with reference [3]. 
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y y 

FIG. 4. (a) Evolution of a single finger in a channel, plotted 
at time intervals of 750, for a = 0.035. (b) and (c) Numeri- 
cal results (lines) and Saffman- Taylor solution (symbols) are 
presented for two values of a that lead to two different finger 
widths (b) A = 0.61 and (c) A = 0.53 (c). 



We have measured the finger width and the finger-tip 
velocity v for different values of the parameter a. The re- 
sults, shown in Fig. 5, are in qualitative agreement with 
experimental results of Pitts [21] and Saffman- Taylor [2] 
and with the numerical results of McLean and Saffman 
[3]. We observe that A tends to one half the channel 
width as the velocity increases [22] . 
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FIG. 5. Finger-width A versus velocity v. Solid line is a 
guide to the eye. Inset shows the dependence of v on the 
parameter a. 



CONCLUSIONS 



A one-sided phase-field model to describe the dynamic 
evolution of unstable interfaces for Hele-Shaw flows in the 
high-viscosity-contrast regime has been proposed. The 
mesoscopic model contains a field equation for a con- 
served order parameter (Model B of Ginzburg-Landau 
phenomenology) and a boundary condition that drives 
the interface out of equilibrium. An asymptotic expan- 
sion to derive the macroscopic equations has been per- 
formed. The phase-field model has been numerically in- 
tegrated and we have analysed different stages of the dy- 
namics. We observe how from a random perturbation 
to the interface, fingers develop. Modes grow and com- 
pete dynamically and the competition ends in a single 
steady-state finger. The width of this finger goes to one 
half of the channel width as the velocity increases. This 
is in agreement with experiments and the existent the- 
ory. We have verify that the shape of the finger tip is in 
good agreement with the parametric solution of Saffman 
and Taylor when the finger width is close to one half of 
the channel width. Also for larger width the shape is 
in qualitative agreement with the fingers found by Mc 
Lean and Saffman. We believe that our model could be 
a useful tool to study situations that cannot be easily 
tackled with traditional methods, like integro-differential 
equations, such as the effect introduced by quenched dis- 
order. 
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